﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Security.Cryptography.X509Certificates;
using System.Text;
using System.Threading.Tasks;

namespace preprocess
{
    internal class Fit
    {
        public List<double> factor = new List<double>(); ///<拟合后的方程系数
        public List<double> Quafactor = new List<double>();
        
        double ssr;                 ///<回归平方和
		double sse;                 ///<(剩余平方和)
		double rmse;                ///<RMSE均方根误差
		List<double> fitedYs = new List<double>();///<存放拟合后的y值，在拟合时可设置为不保存节省内存

        //bool LinearSelect;
        //bool QuaSelect;

        public Fit()
        {
            ssr = 0;
            sse = 0;
            rmse = 0;
            factor.Add(0);
            factor.Add(0);
            for (int i = 0;  i< 3; i++) { Quafactor.Add(0); }
        }
        ~Fit() { }
        /*-------------------------Fit Method 1: Linear Fit------------------------------*/
        /// \brief 直线拟合-一元回归,拟合的结果可以使用getFactor获取，或者使用getSlope获取斜率，getIntercept获取截距
        /// \param x 观察值的x
        /// \param y 观察值的y
        /// \param length x,y数组的长度
        /// \param isSaveFitYs 拟合后的数据是否保存，默认否
        ///
        public bool linearFit(List<double> x, List<double> y, bool isSaveFitYs = false)
        {
            return linearFit(x, y, x.Count, isSaveFitYs);
        }
        public bool linearFit(List<double> x, List<double> y, int length, bool isSaveFitYs = false)
        {
            factor[0] = 0; factor[1] = 0;
            double t1 = 0, t2 = 0, t3 = 0, t4 = 0;
            for (int i = 0; i < length; ++i)
            {
                //Console.WriteLine("data:" + x[i] + "    " + y[i]);
                t1 += x[i] * x[i];
                t2 += x[i];
                t3 += x[i] * y[i];
                t4 += y[i];
            }
            //Console.WriteLine("t:"+t1+"  " + t2 + "  " + t3 + "  " + t4);
            double mt1 = t1 / (length);
            double mt2 = t2 / (length);
            double mt3 = t3 / (length);
            double mt4 = t4 / (length);
            //Console.WriteLine("m:" + mt1 + "  " + mt2 + "  " + mt3 + "  " + mt4);
            factor[0] = (mt3 - mt2 * mt4) / (mt1 - Math.Pow(mt2, 2));
            factor[1] = mt4 - factor[0] * mt2;
            //Console.WriteLine("factor[0]:" + factor[0]);
            //计算误差
            calcError(x, y, length, isSaveFitYs);
            return true;
        }

        /*-------------------------Fit Method 2: Quadratic Fit------------------------------*/
        public bool QuaFit(List<double> x, List<double> y, bool isSaveFitYs = false)
        {
            return QuaFit(x, y, x.Count, isSaveFitYs);
        }

        // Assumed Formular: y = a * x^2 + b * x + c
        // Instruction: Using Least Squares Method (LSM) to get the value of a, b, c
        // 最小二乘法拟合
        public bool QuaFit(List<double> x, List<double> y, int length, bool isSaveFitYs = false)
        {
            //Quafactor[0] = 0; 
            //Quafactor[1] = 0; 
            //Quafactor[2] = 0; // a, b, c 
            // m: mean value; xn : x^n
            // QTV: Quadratic total values: mx_total, my_total, mx2_total, mx3_total, mx4_toatal, mxy_total, mx2y_total
            // QMV: Quadratic mean values: mx, my, mx2, mx3, mx4, mxy, mx2y
            List<double> QTV = new List<double>();
            List<double> QMV = new List<double>();
            for (int i = 0;  i< 7; i++) { QTV.Add(0); }
            for (int i = 0; i< 7; i++) { QMV.Add(0); }


            for (int i = 0; i < length; i++)
            {
                QTV[0] += x[i];
                QTV[1] += y[i];
                QTV[2] += x[i] * x[i];
                QTV[3] += x[i] * x[i] * x[i];
                QTV[4] += x[i] * x[i] * x[i] * x[i];
                QTV[5] += x[i] * y[i];
                QTV[6] += x[i] * x[i] * y[i];
            }
            
            for (int i = 0; i < 7; i++)
            {
                QMV[i] = QTV[i] / length;
            }

            Quafactor[0] = ((QMV[2] - QMV[0] * QMV[1]) * (QMV[3] - QMV[0] * QMV[2]) - ((QMV[6] - QMV[2] * QMV[1]) * (QMV[2] - Math.Pow(QMV[0], 2))))
                        / (Math.Pow((QMV[3] - QMV[0] * QMV[2]), 2) - (QMV[4] - Math.Pow(QMV[2], 2) * QMV[2]) * (QMV[2] - Math.Pow(QMV[0], 2)));
            Quafactor[1] = (QMV[5] - QMV[0] * QMV[1] - Quafactor[0] * (QMV[3] - QMV[0] * QMV[2]))
                        / (QMV[2] - Math.Pow(QMV[0], 2));
            Quafactor[2] = QMV[1] - Quafactor[0] * QMV[2] - Quafactor[1] * QMV[0];

            //////////////////////////////////////////////////////////////////////////
            //计算误差
            //calcQuaError(x, y, length, isSaveFitYs);
            calcQuaError(Quafactor[0], Quafactor[1], x, y, length, isSaveFitYs);
            return true;
        }


        /*-------------------------Fit Method 3: Polynonimal Fit------------------------------*/
        ///
        /// \brief 多项式拟合，拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n
        /// \param x 观察值的x
        /// \param y 观察值的y
        /// \param poly_n 期望拟合的阶数，若poly_n=2，则y=a0+a1*x+a2*x^2
        /// \param isSaveFitYs 拟合后的数据是否保存，默认是
        /// 
        void polyfit(List<double> x, List<double> y, int poly_n, bool isSaveFitYs = true)
        {
            polyfit(x, y, x.Count, poly_n, isSaveFitYs);
        }

        void polyfit(List<double> x, List<double> y, int length, int poly_n, bool isSaveFitYs = true)
        {
            // initial factor list
            factor = new List<double>();
            for (int n = 0; n < poly_n + 1; n++)
            {
                factor.Add(0);
            }


            int i, j;
            //double *tempx,*tempy,*sumxx,*sumxy,*ata;
            List<double> tempx = new List<double>();
            for (int n = 0; n < length; n++)
                tempx.Add(1.0);


            List<double> tempy = new List<double>();
            for (int n = 0; n < y.Count; n++) //??不确定 QF
                tempy.Add(y.Count + length); //??不确定 QF

            List<double> sumxx = new List<double>();
            for (int n = 0; n < poly_n * 2 + 1; n++)
                sumxx.Add(0);

            List<double> ata = new List<double>();
            for (int n = 0; n < (poly_n + 1) * (poly_n + 1); n++)
                sumxx.Add(0);

            List<double> sumxy = new List<double>();
            for (int n = 0; n < (poly_n + 1); n++)
                sumxx.Add(0);

            for (i = 0; i < 2 * poly_n + 1; i++)
            {
                for (sumxx[i] = 0, j = 0; j < length; j++)
                {
                    sumxx[i] += tempx[j];
                    tempx[j] *= x[j];
                }
            }
            for (i = 0; i < poly_n + 1; i++)
            {
                for (sumxy[i] = 0, j = 0; j < length; j++)
                {
                    sumxy[i] += tempy[j];
                    tempy[j] *= x[j];
                }
            }
            for (i = 0; i < poly_n + 1; i++)
                for (j = 0; j < poly_n + 1; j++)
                    ata[i * (poly_n + 1) + j] = sumxx[i + j];
            gauss_solve(poly_n + 1, ata, factor, sumxy);
            //计算拟合后的数据并计算误差
            fitedYs = new List<double>();
            for (int n = 0; n < length; n++)
                fitedYs.Add(0);
            calcError(x, y, length, isSaveFitYs);

        }
        /// 
        /// \brief 获取系数
        /// \param 存放系数的数组
        ///



        public List<double> getFactor() { return factor; }
        public List<double> getQuaFactor() { return Quafactor; }
        /// 
        /// \brief 获取拟合方程对应的y值，前提是拟合时设置isSaveFitYs为true
        ///
        public List<double> getFitedYs() { return fitedYs; }

        /// 
        /// \brief 根据x获取拟合方程的y值
        /// \return 返回x对应的y值
        ///
        public double getY(double x)
        {
            double ans = 0;
            for (int i = 0; i < factor.Count; ++i)
            {
                ans += factor[i] * Math.Pow((double)x, (int)i);
            }
            return ans;
        }

        /// 
        /// \brief 获取斜率
        /// \return 斜率值
        ///
        public double getSlope() { return factor[0]; }
        /// 
        /// \brief 获取截距
        /// \return 截距值
        ///
        public double getIntercept() { return factor[0]; }
        /// 
        /// \brief 剩余平方和
        /// \return 剩余平方和
        ///
        public double getSSE() { return sse; }
        /// 
        /// \brief 回归平方和
        /// \return 回归平方和
        ///
        public double getSSR() { return ssr; }
        /// 
        /// \brief 均方根误差
        /// \return 均方根误差
        ///
        public double getRMSE() { return rmse; }
        /// 
        /// \brief 确定系数，系数是0~1之间的数，是数理上判定拟合优度的一个量
        /// \return 确定系数
        ///
        public double getR_square() { return 1 - (sse / (ssr + sse)); }
        /// 
        /// \brief 获取两个vector的安全size
        /// \return 最小的一个长度
        ///
        public int getSeriesLength(List<double> x, List<double> y)
        {
            return (x.Count > y.Count ? y.Count : x.Count);
        }
        /// 
        /// \brief 计算均值
        /// \return 均值
        ///

        public double Mean(List<double> v)
        {
            return Mean(v, v.Count);
        }

        public double Mean(List<double> v, int length)
        {
            double total = 0;
            for (int i = 0; i < length; ++i)
            {
                total += v[i];
            }
            return (total / length);
        }
        /// 
        /// \brief 获取拟合方程系数的个数
        /// \return 拟合方程系数的个数
        ///
        public int getFactorSize() { return factor.Count; }
        /// 
        /// \brief 根据阶次获取拟合方程的系数，
        /// 如getFactor(2),就是获取y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n中a2的值
        /// \return 拟合方程的系数
        ///
        public double getFactor(int i) { return factor[i]; }
        public double getQuaFactor(int i) { return Quafactor[i];}

        //public double getQuaDydx_start() {return QuaDydx_start;}

        private void calcError(List<double> x, List<double> y, int length, bool isSaveFitYs = true)
        {
            double mean_y = Mean(y, length);
            double yi = 0;
            fitedYs = new List<double>();
            for (int n = 0; n < length; n++)
                fitedYs.Add(0);
            for (int i = 0; i < length; ++i)
            {
                yi = getY(x[i]);
                ssr += ((yi - mean_y) * (yi - mean_y));//计算回归平方和
                sse += ((yi - y[i]) * (yi - y[i]));//残差平方和
                if (isSaveFitYs)
                {
                    fitedYs.Add(yi);
                }
            }
            rmse = Math.Sqrt(sse / length);
        }

        private void calcQuaError(double a, double b, List<double> x, List<double> y, int length, bool isSaveFitYs = true)
        {
            double mean_y = 2 * a * x[(x.Count-1)/ 2] + b; 
            double yi = 0;
            fitedYs = new List<double>();
            for (int n = 0; n < length; n++)
                fitedYs.Add(0);
            for (int i = 0; i < length; ++i)
            {
                yi = getY(x[i]);
                ssr += ((yi - mean_y) * (yi - mean_y));//计算回归平方和
                sse += ((yi - y[i]) * (yi - y[i]));//残差平方和
                if (isSaveFitYs)
                {
                    fitedYs.Add(yi);
                }
            }
            rmse = Math.Sqrt(sse / length);
        }

        private void gauss_solve(int n, List<double> A, List<double> x, List<double> b)
        {
            int i, j, k, r;
            double max;
            for (k = 0; k < n - 1; k++)
            {
                max = Math.Abs(A[k * n + k]); /*find maxmum*/
                r = k;
                for (i = k + 1; i < n - 1; i++)
                {
                    if (max < Math.Abs(A[i * n + i]))
                    {
                        max = Math.Abs(A[i * n + i]);
                        r = i;
                    }
                }
                if (r != k)
                {
                    for (i = 0; i < n; i++)         /*change array:A[k]&A[r] */
                    {
                        max = A[k * n + i];
                        A[k * n + i] = A[r * n + i];
                        A[r * n + i] = max;
                    }
                }
                max = b[k];                    /*change array:b[k]&b[r]     */
                b[k] = b[r];
                b[r] = max;
                for (i = k + 1; i < n; i++)
                {
                    for (j = k + 1; j < n; j++)
                        A[i * n + j] -= A[i * n + k] * A[k * n + j] / A[k * n + k];
                    b[i] -= A[i * n + k] * b[k] / A[k * n + k];
                }
            }

            for (i = n - 1; i >= 0; x[i] /= A[i * n + i], i--)
                for (j = i + 1, x[i] = b[i]; j < n; j++)
                    x[i] -= A[i * n + j] * x[j];
        }

    }
}

